module mini
use params

contains

SUBROUTINE tauchen(n,rho,mu_z,var_e,Z,Pi,pow,m)
!z_{t+1}=(1-rho)*mu_z+rho*z_t+eps_t
!z_{t+1}-mu_z=rho*(z_t-mu_z)+eps_t
!y_{t+1}=rho*y_t+eps_t
!rho: persistence of AR(1) process
!mu_z: Mean of z
!var_e: Variance of epsilon
!n: number of finite states
INTEGER, INTENT(IN) :: n
double precision, intent(in) :: m
DOUBLE PRECISION, INTENT(IN) :: rho,mu_z,var_e, pow
DOUBLE PRECISION, INTENT(OUT) :: Z(n,1),Pi(n,n)
DOUBLE PRECISION :: Y(n,1),sig_z,w
INTEGER :: i,j
DOUBLE PRECISION :: x,mu,sd,pr
mu=0.
sd=SQRT(var_e)

sig_z=SQRT(var_e/(1.-(rho**2.)))

	DO i=1,n	
	Y(i,1)=-m*sig_z+2.*m*sig_z*((i-1.)/(n-1.))**pow
	Z(i,1)=Y(i,1)+mu_z	
	ENDDO
	
	DO i=1,n

	w=Y(2,1)-Y(1,1)
	x=Y(1,1)+w/2.-rho*Y(i,1)
	CALL Ncdf(x,mu,sd,pr)
	Pi(i,1)=pr

	w=Y(2,1)-Y(1,1)	
	x=Y(1,1)-w/2.-rho*Y(i,1)
	CALL Ncdf(x,mu,sd,pr)
	Pi(i,1)=Pi(i,1)-pr

        w=Y(n,1)-Y(n-1,1)	
	x=Y(n,1)+w/2.-rho*Y(i,1)
	CALL Ncdf(x,mu,sd,pr)
	Pi(i,n)=pr
	
        w=Y(n,1)-Y(n-1,1)	
	x=Y(n,1)-w/2.-rho*Y(i,1)
	CALL Ncdf(x,mu,sd,pr)
	Pi(i,n)=Pi(i,n)-pr

		DO j=2,n-1

		w=Y(j+1,1)-Y(j,1)
		x=Y(j,1)+w/2.-rho*Y(i,1)
		CALL Ncdf(x,mu,sd,pr)
		Pi(i,j)=pr

		w=Y(j,1)-Y(j-1,1)		
		x=Y(j,1)-w/2.-rho*Y(i,1)
		CALL Ncdf(x,mu,sd,pr)
		Pi(i,j)=Pi(i,j)-pr
		
		ENDDO
		
	ENDDO
	
	do i=1,n
	Pi(i,:)=Pi(i,:)/sum(Pi(i,:))
	end do

END SUBROUTINE tauchen

SUBROUTINE Ncdf(x,mu,sd,pr)
IMPLICIT NONE
DOUBLE PRECISION, INTENT(IN) :: x,mu,sd
DOUBLE PRECISION, INTENT(OUT) :: pr
DOUBLE PRECISION :: z
DOUBLE PRECISION, PARAMETER :: PI=3.141592653589793238462643383279502884197

z=(x-mu)/sd

	IF (z .GT. 0.) THEN
	pr=.5+.5*SQRT(1.-EXP(-SQRT(PI/8.)*(z**2.)))
	ELSE
	pr=.5-.5*SQRT(1.-EXP(-SQRT(PI/8.)*(z**2.)))
	ENDIF
END SUBROUTINE Ncdf

SUBROUTINE rouwenhorst(n,rho,mu_z,var_e,Z,Pi)
IMPLICIT NONE
!z_{t+1}=(1-rho)*mu_z+rho*z_t+eps_t
!z_{t+1}-mu_z=rho*(z_t-mu_z)+eps_t
!y_{t+1}=rho*y_t+eps_t
!rho: persistence of AR(1) process
!mu_z: Mean of z
!var_e: Variance of epsilon
!n: number of finite states
double precision, INTENT(IN) :: rho,mu_z,var_e
INTEGER, INTENT(IN) :: n
double precision, INTENT(OUT):: Z(n,1),Pi(n,n)
double precision :: p
double precision :: sig_z, psi
double precision :: Y(n,1)
INTEGER :: i,j,k
double precision :: T(n,n,n),T1(n,n,n),T2(n,n,n),T3(n,n,n),T4(n,n,n)
sig_z=SQRT(var_e/(1.-(rho**2.)))
psi=(SQRT(n-1.))*sig_z
p=(1.+rho)/2.

	DO i=1,n
	Y(i,1)=-psi+2.*psi*(i-1.)/(n-1.)
	Z(i,1)=Y(i,1)+mu_z
	ENDDO
	
	if ( n .eq. 1 ) then
	pi(1,1)=1
	z(1,1)=mu_z
	end if

	IF (n.EQ.2) THEN
	Pi(1,1)=p
	Pi(1,2)=1.-p
	Pi(2,1)=1.-p
	Pi(2,2)=p
	ENDIF
T=0.
T1=0.
T2=0.
T3=0.
T4=0.
	IF (n.GT.2) THEN
	T(2,1,1)=p
	T(2,1,2)=1.-p
	T(2,2,1)=1.-p
	T(2,2,2)=p
		DO i=3,n
			DO j=1,i-1
				DO k=1,i-1
				T1(i,j,k)=p*T(i-1,j,k)
				T2(i,j,k+1)=(1.-p)*T(i-1,j,k)
				T3(i,j+1,k)=(1.-p)*T(i-1,j,k)
				T4(i,j+1,k+1)=p*T(i-1,j,k)
				ENDDO		
			ENDDO
		T(i,:,:)=T1(i,:,:)+T2(i,:,:)+T3(i,:,:)+T4(i,:,:)
			DO j=2,i-1
			T(i,j,:)=T(i,j,:)/2.	
			ENDDO
		ENDDO
	
		Pi(:,:)=T(n,:,:)

	ENDIF


END SUBROUTINE rouwenhorst

subroutine sort(x,n)
implicit none
integer, intent(in) :: n
double precision, dimension(n,2), intent(inout) :: x

integer :: i, location
double precision :: y(n)

y(:)=x(:,1)

	do i=1,n-1
	location=findmin(y,i,n)
		call swap(x(i,:),x(location,:))
	y(:)=x(:,1)
	end do

end subroutine sort

function findmin (x,beg,en)
implicit none
double precision, dimension(:), intent(in) :: x
integer, intent(in) :: beg, en

double precision :: minimum
integer :: location, i
integer :: findmin

minimum=x(beg)
location=beg

	do i=beg+1,en
		if ( x(i) .le. minimum ) then
		minimum=x(i)
		location=i
		end if
	end do

findmin=location

end function findmin

subroutine swap(a,b)
implicit none
double precision, intent(inout) :: a(2), b(2)
double precision :: temp(2)
temp=a
a=b
b=temp
end subroutine swap

subroutine percentile (p1,x,nrows,perc)
implicit none
double precision, intent(in) :: p1
integer, intent(in) :: nrows
double precision, intent(in) :: x(nrows,2)
double precision, intent(out) :: perc

integer :: check
double precision :: temp
integer :: i
 
i=1
    check=0
temp=0
    
    do while ( i .le. nrows .and. check .eq. 0 )

    temp=temp+x(i,2)
        
        if ( temp .ge. p1 .and. check .eq. 0 ) then
        check=1
        perc=x(i,1)
        end if
            
    i=i+1
	
    end do
	
    if ( check .eq. 0 .and. rank .eq. 0 ) then
    print*,'did not find p1 in percentile mini.f90 '
        call exit
    end if
	
end subroutine percentile


subroutine onedlin (x,y,n,x1,y1)
implicit none
integer, intent(in) :: n
double precision, intent(in) :: x(n), y(n), x1
double precision, intent(out) :: y1
integer :: point, find
    
    if ( x1 .lt. x(1) ) then
    point=1
    y1=y(point)+(x1-x(point))*(y(point+1)-y(point))/(x(point+1)-x(point))
    else if ( x1 .gt. x(n) ) then
    point=n-1
    y1=y(point)+(x1-x(point))*((y(point+1)-y(point))/(x(point+1)-x(point)))
    
    else
    
    point=1
    find=0
    
        do while ( point .lt. n )

            if ( x1 .ge. x(point) .and. x1 .lt. x(point+1) ) then
            find=1
            y1=y(point)+(x1-x(point))*(y(point+1)-y(point))/(x(point+1)-x(point))
            point=point+n+100		
            end if

        point=point+1
        
        end do

        if ( find .eq. 0 ) then
        print*,'find .eq. 0, onedlin, mini.f90'
            call exit
        end if
        
    end if

    if ( isnan(y1) ) then
    print*,'isnan(y1), onedline, mini.f90'
        call exit
    end if
    
end subroutine onedlin

SUBROUTINE shft(a,b,c,d)
DOUBLE PRECISION , INTENT(OUT) :: a
DOUBLE PRECISION , INTENT(INOUT) :: b,c
DOUBLE PRECISION , INTENT(IN) :: d
	a=b
	b=c
	c=d
END SUBROUTINE shft

function func_sum(x)
implicit none
double precision, intent(in), dimension(4) :: x
double precision :: func_sum

var_z=x(1)
var_eta=x(2)
var_nu=x(3)
rho_eta=x(4)

    call gridearnings

func_sum=(std_log_earnings_change_1_2019-std_log_earnings_change_1)**2+&
(std_log_earnings_change_5_2019-std_log_earnings_change_5)**2+&
(std_log_earnings_perm_2019-std_log_earnings_perm)**2+&
(correlation_lag_1_2019-correlation_lag_t(1))**2+&
(correlation_lag_2_2019-correlation_lag_t(2))**2+&
(correlation_lag_3_2019-correlation_lag_t(3))**2+&
(correlation_lag_4_2019-correlation_lag_t(4))**2+&
(correlation_lag_5_2019-correlation_lag_t(5))**2

end function func_sum

subroutine amoeba (p,y,n,tol,iter_amoeba,func)
implicit none
INTEGER, INTENT(IN) :: n	!dimension of domain for f
INTEGER, INTENT(OUT) :: iter_amoeba
DOUBlE PRECISION, INTENT(OUT), DIMENSION(n) :: p
DOUBLE PRECISION, INTENT(OUT) :: y
DOUBLE PRECISION, INTENT(IN) :: tol
DOUBLE PRECISION, DIMENSION(n+1,n) :: x	
DOUBLE PRECISION, DIMENSION(n) :: x0, x_r, x_e, x_c, x_i	!centroid, reflection, expansion, contraction, reduction
DOUBLE PRECISION, DIMENSION(n+1) :: f 
!DOUBLE PRECISION, PARAMETER :: alpha=1., gamma=2, rho=.5, sigma=scale*.5
DOUBLE PRECISION, PARAMETER :: alpha=1., gamma=2, rho=.5, sigma=.5
!gamma>1    !rho \in [0,.5]
!centroid, reflection, expansion, contraction, reduction  alpha>0, gamma>1, -.5<=rho<0
!recommended alpha=1, gamma=2, rho=-.5, and sigma=.5
INTEGER :: i
INTEGER, PARAMETER :: itmax_amoeba=100
DOUBLE PRECISION :: d

double precision :: func_x_r, func_x_e, func_x_c

    interface 
        function func (y)
        double precision, intent(in), dimension(4) :: y
        double precision :: func
        end function func
    end interface 

d=1.
iter_amoeba=0

!set simplex here
    do iter=1,n+1
    x(iter,:)=rn(iter_calibrate,rank+1,iter,:)
    f(iter)=func(x(iter,:))
    end do

    call step_1(x,f,n)	!ascending order

    do while ( iter_amoeba .lt. itmax_amoeba .and. d .gt. tol )

        call step_2(x0,x,n) !centroid
        
    x_r(:)=x0(:)+alpha*(x0(:)-x(n+1,:)) !alpha > 0
    
        do i=1,n
        x_r(i)=max(1e-6,x_r(i))
        end do

        do i=1,n
        x_r(i)=min(one*.999,x_r(i))
        end do
        
    func_x_r=func(x_r)
    
!         if ( rank .eq. 0 ) then
!         print*,'reflection',iter_amoeba,d,f(1)
!         end if
    
        if ( f(1) .le. func_x_r .and. func_x_r .lt. f(n) ) then

        !reflection
        x(n+1,:)=x_r(:)
        f(n+1)=func_x_r
            
        else if ( func_x_r .lt. f(1) ) then	

        !expansion
!         x_e(:)=x0(:)+gamma*alpha*(x0(:)-x(n+1,:))
        x_e(:)=x0(:)+gamma*(x_r(:)-x0(:))

            do i=1,n
            x_e(i)=max(1e-6,x_e(i))
            end do

            do i=1,n
            x_e(i)=min(one*.999,x_e(i))
            end do        
        
        func_x_e=func(x_e)
!         iter_amoeba=iter_amoeba+1
        
!             if ( rank .eq. 0 ) then
!             print*,'expansion',iter_amoeba,d,f(1)
!             end if
		
            if ( func_x_e .lt. func_x_r ) then
            x(n+1,:)=x_e(:)
            f(n+1)=func_x_e
            else if ( func_x_e .ge. func_x_r ) then
            x(n+1,:)=x_r(:)
            f(n+1)=func_x_r
            end if
            
        else if ( func_x_r .ge. f(n) ) then

        !contraction	
        
            if ( func_x_r .lt. f(n+1) ) then
            
            x_c(:)=x0(:)+rho*(x_r(:)-x0(:))   !0 < rho <= 0.5

                do i=1,n
                x_c(i)=max(1e-6,x_c(i))
                end do

                do i=1,n
                x_c(i)=min(one*.999,x_c(i))
                end do        
            
            func_x_c=func(x_c)            
!             iter_amoeba=iter_amoeba+1

!             if ( rank .eq. 0 ) then
!             print*,'contraction',iter_amoeba,d,f(1)
!             end if
        
                if ( func_x_c .lt. func_x_r ) then

                x(n+1,:)=x_c(:)
                f(n+1)=func_x_c

                else if ( func_x_c .ge. func_x_r ) then
                
                !reduction            
                    do i=2,n+1                
                    
                    x(i,:)=x(1,:)+sigma*(x(i,:)-x(1,:))     
                    
                    x(i,:)=max(1e-6,x(i,:))
                    x(i,:)=min(one*.999,x(i,:))
                
                    f(i)=func(x(i,:))
!                     iter_amoeba=iter_amoeba+1
                
!                     if ( rank .eq. 0 ) then
!                     print*,'reduction',iter_amoeba,d,f(1)
!                     end if
                
                    end do                
                
                end if

            else if ( func_x_r .ge. f(n+1) ) then

            x_c(:)=x0(:)+rho*(x(n+1,:)-x0(:))   !0 < rho <= 0.5
            
                do i=1,n
                x_c(i)=max(1e-6,x_c(i))
                end do

                do i=1,n
                x_c(i)=min(one*.999,x_c(i))
                end do        
            
            func_x_c=func(x_c)                                
!             iter_amoeba=iter_amoeba+1

                if ( func_x_c .lt. f(n+1) ) then

                x(n+1,:)=x_c(:)
                f(n+1)=func_x_c
            
                else if ( func_x_c .ge. f(n+1) ) then

                !reduction            
                    do i=2,n+1                
                    
                    x(i,:)=x(1,:)+sigma*(x(i,:)-x(1,:))     

                    x(i,:)=max(1e-6,x(i,:))
                    x(i,:)=min(one*.999,x(i,:))
                    
                    f(i)=func(x(i,:))
!                     iter_amoeba=iter_amoeba+1
                
!                     if ( rank .eq. 0 ) then
!                     print*,'reduction',iter_amoeba,d,f(1)
!                     end if
                
                    end do

                end if
                
            end if

        end if

        call step_1(x,f,n)	!ascending order

    d=(maxval(x(:,1))-minval(x(:,1)))/minval(x(:,1))
    iter_amoeba=iter_amoeba+1
   
        do i=2,n
        
            if ( (maxval(x(:,i))-minval(x(:,i)))/minval(x(:,i)) .ge. d ) then
            d=(maxval(x(:,i))-minval(x(:,i)))/minval(x(:,i))
            end if
            
        end do

    end do

    if ( minloc(f,1) .ne. 1 ) then
    print*,'mistake in amoeba',f
        call exit
    end if
    
p(:)=x(1,:)
y=f(1)

!    if ( iter_amoeba .ge. itmax_amoeba ) then
!    print*,'iterations exceeded at amoeba',d
!    end if

end subroutine amoeba

SUBROUTINE step_2 (x0,x,n)
!centroid
IMPLICIT NONE
INTEGER, INTENT(IN) :: n
DOUBLE PRECISION, INTENT(IN), DIMENSION(n+1,n) :: x
DOUBLE PRECISION, INTENT(INOUT), DIMENSION(n) :: x0
INTEGER :: i, j

x0=0

    DO i=1,n
        DO j=1,n	
        x0(i)=x0(i)+x(j,i)
        ENDDO
    ENDDO

x0=x0/n

END SUBROUTINE step_2

SUBROUTINE step_1 (x,f,n)
!ascending order
IMPLICIT NONE
INTEGER, INTENT(IN) :: n
DOUBLE PRECISION, INTENT(INOUT), DIMENSION(n+1,n) :: x
DOUBLE PRECISION, INTENT(INOUT), DIMENSION(n+1) :: f
INTEGER :: i, j, loc
DOUBLE PRECISION :: swap_1, swap_2(n)
double precision :: minimum

    DO i=1,n+1
        IF ( ISNAN(f(i)) ) THEN
        PRINT*,'Nan at step_1'
            CALL EXIT
        ENDIF
    ENDDO

    DO i=1,n
    
    loc=i
    minimum=f(loc)
    
        DO j=i+1,n+1
        
            IF ( f(j) .LT. minimum ) THEN
            loc=j
            minimum=f(j)
            ENDIF
            
        ENDDO
		
    swap_1=f(i)
    swap_2(:)=x(i,:)
    
    f(i)=f(loc)
    x(i,:)=x(loc,:)
    
    f(loc)=swap_1
    x(loc,:)=swap_2(:)
    
    ENDDO
    
END SUBROUTINE step_1

end module mini 
